NASA/CR-2002-211648 
ICASE  Report  No.  2002-13 


Analysis  of  Boundary  Conditions  for  Factorizable 
Discretizations  of  the  Euler  Equations 


Boris  Diskin 

ICASE,  Hampton,  Virginia 
James  L.  Thomas 

NASA  Langley  Research  Center,  Hampton,  Virginia 


May  2002 


Report  Documentation  Page 


Report  Date 

Report  Type 

Dates  Covered  (from...  to) 

00MAY2002 

N/A 

- 

Title  and  Subtitle 

Contract  Number 

Analysis  ot  Boundary  Conditions  tor  Factorizable 
Discretizations  of  the  Euler  Equations 

Grant  Number 

Program  Element  Number 

Author(s) 

Project  Number 

Boris  Diskin, James  L.  Thomas 

Task  Number 

Work  Unit  Number 

Performing  Organization  Name(s)  and  Address(es) 

ICASE  Mail  Stop  132C  NASA  Langley  Research  Center 
Hampton,  VA  23681-2199 

Performing  Organization  Report  Number 

Sponsoring/Monitoring  Agency  Name(s)  and 

Sponsor/Monitor’s  Acronym(s) 

Address(es) 

Sponsor/Monitor’s  Report  Number(s) 

Distribution/Availability  Statement 

Approved  for  public  release,  distribution  unlimited 

Supplementary  Notes 

ICASE  Report  No.  2002-13 

Abstract 

see  report 

Subject  Terms 

Report  Classification 

unclassified 

Classification  of  this  page 

unclassified 

Classification  of  Abstract 

unclassified 

Limitation  of  Abstract 

SAR 

Number  of  Pages 

31 

The  NASA  STI  Program  Office  ...  in  Profile 


Since  its  founding,  NASA  has  been  dedicated 
to  the  advancement  of  aeronautics  and  space 
science.  The  NASA  Scientific  and  Technical 
Information  (STI)  Program  Office  plays  a  key 
part  in  helping  NASA  maintain  this 
important  role. 

The  NASA  STI  Program  Office  is  operated  by 
Langley  Research  Center,  the  lead  center  for 
NASA’s  scientific  and  technical  information. 
The  NASA  STI  Program  Office  provides 
access  to  the  NASA  STI  Database,  the 
largest  collection  of  aeronautical  and  space 
science  STI  in  the  world.  The  Program  Office 
is  also  NASA’s  institutional  mechanism  for 
disseminating  the  results  of  its  research  and 
development  activities.  These  results  are 
published  by  NASA  in  the  NASA  STI  Report 
Series,  which  includes  the  following  report 
types: 

•  TECHNICAL  PUBLICATION.  Reports  of 
completed  research  or  a  major  significant 
phase  of  research  that  present  the  results 
of  NASA  programs  and  include  extensive 
data  or  theoretical  analysis.  Includes 
compilations  of  significant  scientific  and 
technical  data  and  information  deemed 

to  be  of  continuing  reference  value.  NASA’s 
counterpart  of  peer-reviewed  formal 
professional  papers,  but  having  less 
stringent  limitations  on  manuscript 
length  and  extent  of  graphic 
presentations. 

•  TECHNICAL  MEMORANDUM. 

Scientific  and  technical  findings  that  are 
preliminary  or  of  specialized  interest, 
e.g.,  quick  release  reports,  working 
papers,  and  bibliographies  that  contain 
minimal  annotation.  Does  not  contain 
extensive  analysis. 

•  CONTRACTOR  REPORT.  Scientific  and 
technical  findings  by  NASA-sponsored 
contractors  and  grantees. 


•  CONEERENCEPUBEICATIONS. 
Collected  papers  from  scientific  and 
technical  conferences,  symposia, 
seminars,  or  other  meetings  sponsored  or 
cosponsored  by  NASA. 

•  SPECIAE  PUBEICATION.  Scientific, 
technical,  or  historical  information  from 
NASA  programs,  projects,  and  missions, 
often  concerned  with  subjects  having 
substantial  public  interest. 

•  TECHNICAE  TRANSEATION.  English- 
language  translations  of  foreign  scientific 
and  technical  material  pertinent  to 
NASA’s  mission. 

Specialized  services  that  complement  the 
STI  Program  Office’s  diverse  offerings  include 
creating  custom  thesauri,  building  customized 
data  bases,  organizing  and  publishing 
research  results  .  .  .  even  providing  videos. 

Eor  more  information  about  the  NASA  STI 
Program  Office,  see  the  following: 

•  Access  the  NASA  STI  Program  Home 
Page  at  http://www.sti.nasa.gov 

•  Email  your  question  via  the  Internet  to 
help@sti.nasa.gov 

•  Eax  your  question  to  the  NASA  STI 
Help  Desk  at  (301)  621-0134 

•  Telephone  the  NASA  STI  Help  Desk  at 
(301)  621-0390 

•  Write  to: 

NASA  STI  Help  Desk 

NASA  Center  for  AeroSpace  Information 

7121  Standard  Drive 

Hanover,  MD  21076-1320 


NASA/CR-2002-211648 
ICASE  Report  No.  2002-13 


Analysis  of  Boundary  Conditions  for  Factorizable 
Discretizations  of  the  Euler  Equations 


Boris  Diskin 

ICASE,  Hampton,  Virginia 
James  L.  Thomas 

NASA  Langley  Research  Center,  Hampton,  Virginia 


ICASE 

NASA  Langley  Research  Center 
Hampton,  Virginia 

Operated  by  Universities  Space  Research  Association 


Prepared  for  Langley  Research  Center 
under  Contract  NAS  1-97046 


May  2002 


Available  from  the  following: 

NASA  Center  for  AeroSpace  Information  (CASI) 
7121  Standard  Drive 
Hanover,  MD  21076-1320 
(301)  621-0390 


National  Technical  Information  Service  (NTIS) 
5285  Port  Royal  Road 
Springfield,  VA  22161-2171 
(703)  487-4650 


ANALYSIS  OF  BOUNDARY  CONDITIONS  FOR  FACTORIZABLE  DISCRETIZATIONS 

OF  THE  EULER  EQUATIONS 

BORIS  DISKIN*  AND  JAMES  L.  THOMASt 


Abstract.  In  this  article,  several  sets  of  boundary  conditions  for  factorizable  schemes  corresponding 
to  the  steady-state  compressible  Euler  equations  are  evaluated.  The  analyzed  model  is  a  one-dimensional 
constant-coefficient  problem.  Numerical  tests  have  been  performed  for  a  fully  subsonic  quasi-one-dimensional 
flow  in  a  convergent /divergent  channel. 

This  paper  focuses  on  the  effect  of  boundary-condition  equations  on  stability  and  accuracy  of  the  discrete 
solutions.  Explicit  correspondence  between  solutions  and  boundary  conditions  is  established  through  a 
boundary-condition-sensitivity  (BCS)  matrix.  The  following  new  findings  are  reported: 

(1)  Examples  of  stable  discrete  problems  contradicting  a  wide-spread  belief  that  employment  of  a  one- 
order-lower  approximation  schemes  in  an  0(/i)-small  region  does  not  affect  the  overall  accuracy  order 
of  the  solution  have  been  found  and  explained.  Such  counterexamples  can  only  be  constructed  for 
systems  of  differential  equations.  Eor  scalar  equations,  the  conventional  wisdom  is  correct. 

(2)  A  negative  effect  of  overspecified  (although,  exact)  boundary  conditions  on  accuracy  and  stability 
of  the  solution  has  been  observed  and  explained. 

(3)  Sets  of  practical  boundary  conditions  for  factorizable  schemes  providing  stable  second-order  accurate 
solutions  have  been  formulated.  These  schemes  belong  to  a  family  of  second-order  schemes  requiring 
second-order  accuracy  for  some  numerical-closure  boundary  conditions. 

Key  words,  compressible  Euler  equations,  factorizable  discretization  scheme,  practical  boundary  con¬ 
ditions,  I-stability,  B-stability 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  An  efficiency  goal  for  the  new  generation  of  multigrid  flow  solvers  is  to  arrive  at 
solutions  of  the  governing  system  of  equations  in  a  total  computational  work  that  is  a  small  (less  than  10) 
multiple  of  the  operation  count  in  one  target-grid  residual  evaluation.  Such  efficiency  is  defined  as  textbook 
multigrid  efficiency  (TME)  [4,  5].  TME  has  been  achieved  for  elliptic  problems  long  ago  [1,  2,  3,  6].  The 
difficulties  associated  with  extending  TME  for  solution  of  the  Navier-Stokes  equations  relate  to  the  fact  that 
the  Navier-Stokes  equations  are  a  system  of  coupled  nonlinear  equations  that  is  not  fully  elliptic,  even  for 
subsonic  Mach  numbers,  but  contains  hyperbolic  partitions.  TME  for  the  Navier-Stokes  simulations  can  be 
achieved  if  different  factors  contributing  to  the  system  could  be  separated  and  treated  optimally,  e.g.,  by 
multigrid  for  elliptic  factors  and  by  downstream  marching  for  hyperbolic  factors.  An  efficient  way  to  separate 
the  factors  is  the  distributed  relaxation  approach  proposed  in  [2,  3].  The  general  framework  for  achieving 
TME  in  large-scale  computational  fluid  dynamics  (CED)  applications  has  been  discussed  in  [8,  20]. 

Design  of  a  distributed  relaxation  scheme  for  Navier-Stokes  systems  can  be  significantly  simplified  if  the 
target  discretization  possesses  two  properties: 
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(1)  The  discretization  of  the  corresponding  principally  linearized  system  L  is  factorizable  [3,4,  12,  16,  17] , 
i.e.,  the  discrete  system  determinant  can  be  represented  as  a  product  of  discrete  scalar  factors,  each 
of  them  approximating  a  corresponding  factor  of  the  determinant  of  the  differential  Navier-Stokes 
equations. 

(2)  The  obtained  scalar  factor  discretizations  are  stable,  easily  solvable,  and  reflect  the  physical  anisotropies. 

Some  well-known  discrete  schemes  are  naturally  factorizable,  e.g.,  the  staggered-grid  discretization 

scheme  for  incompressible  flow  dating  back  to  the  mid  60’s  [13,  14,  15].  However,  the  majority  of  dis¬ 
crete  schemes  in  current  use  are  not  factorizable.  Even  for  factorizable  schemes,  the  discretizations  obtained 
for  scalar  factors  of  the  determinant  are  not  often  stable  or/and  have  wrong  strong-anisotropy  directions. 
The  search  for  new  factorizable  discretization  schemes  is  chiefly  motivated  by  the  need  to  derive  discrete 
schemes  with  the  resulting  discretizations  of  scalar  factors  satisfying  some  desired  properties  (e.g.,  stability, 
correct  alignment  with  the  physical  anisotropies,  compactness,  availability  of  an  efficient  relaxation  scheme, 
etc.)  Recently,  two  families  of  factorizable  discretization  schemes  for  the  compressible  Euler  equations  have 
emerged  [12,  16]. 

The  original  paper  [12]  has  demonstrated  TME  in  solving  factorizable  discretizations  of  the  compressible 
Euler  equations  with  solution  values  at  the  boundary  and,  wherever  necessary,  outside  of  the  target  compu¬ 
tational  domain  overspecified  from  the  known  exact  solution  of  the  differential  problem.  This  formulation 
is  very  attractive  because  it  bypasses  difficulties  associated  with  boundary  conditions  and  facilitates  TME 
demonstration.  However,  overspecified  boundary  conditions  are  impractical;  usually,  the  boundary  data  are 
not  available  outside  of  the  computational  domain  and  even  at  the  boundaries  the  data  are  incomplete  and/or 
accurate  only  to  some  finite  precision  that  might  depend  on  the  mesh.  Multigrid  considerations  make  the 
overspecified  boundary  conditions  even  less  viable:  on  coarse  grids,  specification  of  data  located  far  beyond 
the  target  computational  domain  is  required.  This  paper  studies  more  practical  boundary  conditions  for  the 
factorizable  discretization  schemes  introduced  in  [12]. 

In  general,  boundary  conditions  for  a  discrete  problem  may  include  two  components:  the  physical  bound¬ 
ary  conditions  and  the  numerical-closure  equations.  The  former  are  discrete  approximations  to  the  boundary 
conditions  of  the  target  differential  problem.  The  latter  are  special  discrete  approximations  (different  from 
the  interior  approximations)  to  the  differential  equations  usually  defined  within  an  0(/i)-small  boundary 
neighborhood. 

This  paper  focuses  on  the  effect  of  boundary-condition  equations  on  stability  of  the  discrete  solutions. 
Eor  steady-state  discrete  problems,  there  are  two  relevant  types  of  stability:  (A)  stability  of  difference 
approximations  to  differential  operators  and  (B)  stability  of  discrete  solutions  with  respect  to  boundary  data. 
Because  the  first  type  of  stability  describes  the  properties  of  the  interior  approximations  it  is  referenced  in 
this  paper  as  I-stability;  the  second  type  of  stability  is  referenced  as  B-stability. 

Eor  a  discrete  operator  L**  approximating  the  differential  constant-coefficient  operator  L,  the  property 
of  I-stability  can  be  defined  as  follows:  Let  the  symbols  of  the  differential  operator  L[u]  and  the  discrete 
operator  L**]?!]  be  defined  as  L(q:)  =  e“*("'®)L[e*("'®)]  and  L*’(w)  =  respectively,  where 

e  =  (x,y,z)  is  a  coordinate  vector,  a  =  (ax,ay,az)  are  frequencies  of  a  continuous  Eourier  component, 
j  =  (jx,jy,jz)  are  the  grid  indexes,  and  w  =  (ujx,iOy,uJz),0  <  jw^,],  jwj,],  <  tt  are  normalized  Eourier 
frequencies. 

Definition  of  I-stability:  The  discrete  operator  L**  is  an  I-stable  approximation  to  the  differential  operator 
L  if  L*’(w)  =  0  implies  that  L(w//i)  =  0. 

This  definition  of  I-stability  requires  from  the  discrete  operator  to  have  no  global  (Eourier-component) 
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solutions  to  the  homogeneous  discrete  equations  besides  those  that  are  also  solutions  of  the  corresponding 
homogeneous  differential  equations.  Most  higher-than-first-order  difference  approximations  admit  some  spu¬ 
rious  solutions  that  do  not  approximate  the  solutions  of  the  differential  equations.  For  I-stable  difference 
approximations,  these  spurious  solutions  should  be  local,  i.e.,  (nearly)  vanishing  beyond  some  neighborhood 
that  shrinks  to  a  zero-measure  set  as  the  mesh  size  h  tends  to  zero.  In  boundary-value  problems,  the  spurious 
solutions  of  I-stable  difference  approximations  are  usually  exponential  functions  of  the  grid  indexes.  These 
functions  represent  the  discrete  boundary  layers  and  fast  decrease  from  the  boundary  toward  the  interior. 
A  wide  (2/i)  central-differencing  approximation  to  the  convection  operator  {a  ■  Sj),  where  d  is  a  constant 
vector,  represents  an  I-unstable  approximation  because  the  corresponding  homogeneous  discrete  equation 
admits  several  solutions:  a  constant  corresponding  to  a  physical  constant  solution  of  the  differential  equation 
and  functions  oscillating  with  the  highest  normalized  frequency  tt  in  some  (or  all)  coordinate  grid  directions. 
The  latter  functions  represent  global  spurious  discrete  solutions  with  the  unit  amplitude  everywhere  in  the 
interior.  These  discrete  solutions  do  not  approximate  solutions  of  the  differential  equation.  Global  spurious 
solutions  are  not  necessarily  highly  oscillating.  Varying  the  Mach  number  parameter  in  the  Scheme  #  1  from 
[12],  one  can  easily  construct  examples  of  I-unstable  approximations  with  smooth  global  spurious  solutions. 

The  notion  of  I-stability  is  a  generalization  of  the  notion  of  /i-ellipticity  widely  used  in  multigrid  theory. 
A  discrete  operator  is  /i-elliptic  if  the  corresponding  discrete  homogeneous  equation  does  not  admit  highly 
oscillatory  solutions.  If  the  target  differential  solutions  are  smooth  on  a  given  grid,  I-stability  implies  h- 
ellipticity. 

For  constant-coefficient  operators,  a  difference  approximation  is  I-stable  if  the  number  of  characteristic- 
polynomial  roots  with  the  unit  amplitude  is  equal  to  the  number  of  linearly  independent  solutions  of  the 
homogeneous  differential  equation.  The  number  of  roots  with  amplitudes  differing  from  one  may  be  arbitrary. 
The  discrete  solutions  corresponding  to  the  roots  with  the  unit  amplitude  are  referred  as  physical  solutions. 
The  physical  solutions  are  associated  with  the  physical  boundary  conditions.  The  coefficients  of  the  spurious 
solutions  are  strongly  influenced  by  the  numerical-closure  equations.  All  the  difference  approximations 
analyzed  in  this  paper  are  I-stable. 

Let  us  define  the  boundary  data  as  known  quantities  used  in  formulation  boundary  conditions,  e.  g., 
solution  values  at  the  boundary  (Dirichlet  data),  solution  derivatives  at  the  boundary  (Neumann  data),  the 
source  functions  for  the  interior  equations  (numerical-closure  data),  etc. 

Definition  of  B-stability:  A  discrete  scheme  approximating  a  differential  problem  is  stable  with  respect  to 
boundary  data  (B-stable)  if  for  any  fixed  (independent  of  mesh  size  h)  perturbation  of  boundary  data,  the 
sequence  of  discrete  solutions  converges  as  h  tends  to  zero. 

B-stability  implies  that  the  discrete  scheme  remains  meaningful  even  if  the  boundary  data  are  not 
precisely  specified.  Even  I-stable  difference  approximations  with  physically  realistic  boundary  conditions 
may  become  B-unstable.  Moreover,  B-stability  may  depend  on  norms  in  which  solution  convergence  is 
considered.  Examples  of  discrete  schemes  for  boundary-value  problems  that  are  B-stable  in  some  common 
integral  norms  but  B-unstable  in  the  Loo  norm  are  shown  in  Section  6. 

This  paper  also  represents  an  attempt  to  formalize  our  intuitive  understanding  of  practical  discrete 
boundary  conditions.  Assuming  that  the  target  differential  problem  has  a  reasonable  set  of  physical  boundary 
conditions  and  is  well  posed,  the  four  requirements  for  practical  discrete  boundary  conditions  supplementing 
an  I-stable  discrete  approximation  in  the  interior  are  formulated  as  following: 

(1)  Location  requirement:  The  discretization  of  physical  boundary  conditions  should  be  specified  at  the 
boundary. 
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(2)  Well-posedness  requirement:  The  obtained  discrete  problem  should  be  well  posed,  i.e.,  possess  a 
unique  solution  (on  a  given  grid)  continuously  dependent  on  the  problem  data. 

(3)  B-stahility  requirement:  The  discrete  problem  should  be  B-stable. 

(4)  Aeeuraey  requirement:  The  aeeuraey  of  the  discrete  solutions  should  be  determined  by  the  approx¬ 
imation  order  of  the  interior  discrete  equations,  i.e.,  should  not  deteriorate  because  of  boundary 
conditions,  either  the  physical  or  numerical-closure  type. 

The  overspecified  boundary  conditions,  while  often  leading  to  well-posed  discrete  formulations,  obviously 
violate  the  location  requirement  (1);  what  is  less  obvious,  in  some  cases,  they  also  violate  the  B-stability  (3) 
and  accuracy  (4)  requirements.  Examples  are  shown  in  Section  6. 

A  detailed  discrete  analysis  is  employed  in  this  paper  to  evaluate  several  sets  of  boundary  conditions.  The 
model  is  a  one-dimensional  constant-coefficient  problem  corresponding  to  the  compressible  Euler  equations, 
although  the  methodology  applies  to  general  systems.  The  interior  discretizations  are  I-stable  factorizable 
schemes  derived  in  [12].  Explicit  correspondence  between  solutions  and  boundary  conditions  is  established 
through  a  boundary-condition-sensitivity  (BCS)  matrix.  Analysis  of  coefficients  of  the  BCSmatrix  provides 
reliable  predictions  for  B-stability  and  accuracy  of  the  discrete  solutions.  The  following  new  findings  are 
reported: 

(1)  Examples  of  I-stable  discrete  problems  contradicting  a  wide-spread  belief  that  employment  of  a 
one-order-lower  approximation  scheme  as  numerical-closure  equations  in  an  0(/i)-small  region  does 
not  affect  the  overall  accuracy  order  of  the  solution  have  been  found  and  explained.  Such  coun¬ 
terexamples  can  only  be  constructed  for  systems  of  differential  equations.  Eor  scalar  equations,  the 
conventional  wisdom  is  correct. 

(2)  A  negative  effect  of  overspecified  (although,  exact)  boundary  conditions  on  accuracy  and  B-stability 
of  the  solution  has  been  observed  and  explained. 

(3)  Sets  of  practical  boundary  conditions  for  factorizable  schemes  of  [12]  that  provide  B-stable  second- 
order  accurate  solutions  have  been  formulated.  These  schemes  belong  to  a  family  of  second-order 
schemes  requiring  second-order  accuracy  for  some  numerical-closure  boundary  conditions. 

The  material  in  this  paper  has  been  organized  in  the  following  order:  Section  2  formulates  the  model 
differential  constant-coefficient  problem.  A  family  of  factorizable  discretizations  for  the  model  problem  is 
outlined  in  Section  3.  The  discrete  solution  structure  is  analyzed  in  Section  4.  Section  5  introduces  and 
illustrates  the  analysis  of  discrete  boundary  conditions.  Section  6  defines  and  analyzes  several  sets  of  bound¬ 
ary  conditions  employed  for  two  factorizable  discretizations.  The  results  of  numerical  tests  performed  for 
quasi-one-dimensional  compressible  subsonic  flow  in  a  convergent-divergent  channel  are  reported  in  Section  7. 
Section  8  contains  concluding  remarks  as  well  as  directions  for  future  research. 

2.  Model  Problem.  The  set  of  the  quasi-one-dimensional  nonconservative  Euler  equations  is  given  by 

udxU  +  jdxP  =  0, 

(2.1)  pc^dxU  +  udxP  =  -jpu^, 

{-/ -  l)ed:^u  +  ud:^e  = 

where  (t{x)  is  the  area  distribution.  The  pressure  p,  internal  (thermal)  energy  e,  density  p,  and  sound  speed 
c  are  related  as 


P={'1-  l)pe, 


(2.2) 

(2.3) 
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where  7  is  the  ratio  of  specific  heats. 

The  corresponding  model  problem 

/  udx  jdx  0 

(2.4)  pc^dx  udx  0 

\  (7  -  l)edx  0  ud, 


u  '' 
P 

e  / 


h{x) 
h{x)  ) 


assumes  that  the  coefficients  u,p,c,  and  e  are  constants  unrelated  to  the  unknown  functions  (u,p,e).  In 
the  subsonic  regime,  a  natural  set  of  boundary  conditions  is  u  and  e  specified  at  the  inflow  boundary  and  p 
specified  at  the  outflow  boundary.  With  this  set  of  boundary  conditions  the  problem  (2.4)  is  well  posed. 

The  third  equation  is  decoupled  from  the  other  equations.  Thus,  for  the  purpose  of  analysis,  we  can 
focus  on  the  system  of  two  constant-coefficient  equations 


(2.5) 


Lq  =  f, 


where 


(2.6) 


/  udx  ^dx  \ 

y  pc^dx  udx  J  ’ 


q  =  {u,p)'^ ,  and  f  =  (/i,  /2)^  The  determinant  of  the  matrix  operator  L  is  the  full-potential  operator 


(2.7)  =  {u^  -  c^)  dxx- 

3.  Factorizable  Discretizations.  In  this  section,  we  briefly  consider  derivation  of  factorizable  dis¬ 
cretizations  for  the  model  problem  (2.5).  Such  schemes  for  the  three-dimensional  Euler  equations  have  been 
derived  in  [12];  generalization  to  conservative  schemes  has  been  discussed  in  [11]. 

The  starting  point  is  a  “basic”  collocated-grid  discretization,  L^’basic,  for  the  matrix  operator  L  of  (2.6) 
defined  as 


(3.1) 


T  h  _ 

^  basic  — 


ud^  i^c 

p 

pc^d^  ud^ 


where  the  discrete  derivatives,  d^,d^,  and  d'^  are  second-order  accurate  upwind  (upwind-biased),  central, 
and  downwind  (downwind-biased)  difference  approximations,  respectively.  The  determinant  of  L^’basic  is  a 
discrete  approximation  to  the  full-potential  operator  given  by 

(3.2)  u'^d^d'^  -  {d^f  . 


This  discrete  approximation  is  not  I-stable  for  subsonic  Much  numbers  {M  =  u/c  <  1)  and  has  an  incorrect 
(streamwise)  direction  of  strong  coupling  for  near-sonic  Much  numbers  {M  1).  An  I-stable  discrete 
approximation  for  the  full-potential  operator  (2.7)  would  be 


(3.3) 


=  (^2  _  ^2^  Q 


h 

XX  ’ 


where  dxx  is  an  I-stable  approximation  to  the  second  derivative.  A  way  to  achieve  I-stability  for  the  discrete 
full-potential  factor  is  to  replace  the  discretization  with  -I-  .4^,  so  that  the  discrete  full-potential 
operator  is  transformed  to  a  desired  (I-stable)  one.  This  transformation  implies  that 

A’^={ud^y^ 


-  [U^d^d^- 


(3.4) 

(3.5) 


'jyh  _ 
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where  is  a  desired  discretization  of  the  full-potential  operator.  Assuming  that  is  second-order  accurate, 
is  0(/i^)-small,  and  the  overall  second-order  discretization  accuracy  is  not  compromised.  The  operator 
(ud^)  "  is  a  nonlocal  operator  and  its  introduction  can  be  effected  through  a  new  auxiliary  variable  and 
a  new  discrete  equation  =  "D^^-  Finally,  a  family  of  factorizable  I-stable  discretizations  for  (2.6)  is 

formulated  as 

n(q*’)i  :=  ud^u^j  +  =  fj, 

(3.6)  r2(q*’)j  :=  =  0, 

»’3(q*')i  :=  pc^d^Uj  +'tpj  +  ud'^Pj  =  fj, 
where  q*’^  =  (u^,  V'j  ,p^)^,  fj  =  fi(jh),  and  fj  =  f2(jh) 

4.  Solutions.  Let  the  exact  solution  of  the  differential  problem  (2.5)  defined  on  the  interval  x  €  [0, 1] 
be 


(4.1) 


Qexact  — 


‘^exact 

Pexact 


where  a  is  an  arbitrary  frequency.  Then 


Cu 

c„ 


iax 

^  5 


(4.2) 


f  fiix)  \  ^  f  h 
V  h{x)  )  \  f2 


iax 


The  corresponding  system  of  discrete  equations  is 


uCu  +  j;Cp  \  . 

pc^Cu  +  uCp  j 


(4.3) 


ud^ 

0 

\  (  u’^ 

0 

ud^ 

I 

pc^d^ 

1 

V  p'' 

) 

J  3 
0 

J  3 


where  j  =  1, 2, . . . ,  A"  —  1,  N  =  1/h,  fj  =  fi{jh),  and  fj  =  f2{jh).  A  discrete  representation  of  the  exact 
solution  (4.1)  on  a  grid  with  mesh  size  h  is  given  by 


(4-4)  qiact  = 

^  ^exact 
ih^ 

!  V^exact  j 

(  \ 
0 

I 

\  ^^exact  / 

V  Cp  y 

where  to  =  ah  is  a  normalized  frequency.  The  system  (4.3)  is  subject  to  physical  boundary  conditions 


(4.5) 


Uo  =  Cu,  tpo  =  0,  PN  =  CpC°‘. 


Depending  on  the  choice  of  difference  operators,  some  numerical-closure  equations  may  be  required  to  com¬ 
plete  the  discrete  formulation.  The  solution  to  the  problem  (4.3),  (4.5)  can  be  sought  as  a  combination 
of  a  particular  solution  to  the  nonhomogeneous  system  of  equations  (4.3)  and  the  general  solution  to  the 
corresponding  homogeneous  problem  (3.6). 

A  particular  solution  to  (4.3)  can  be  found  in  the  form 


(  ^^par  ^ 

(4-6)  qpar  = 

j  V'par  j 

= 

V  Ppar  / 

K  p  J 

6 


(4.7) 


where 


(4.8) 


L^A)  = 


/  ud^{X)  0 


ia^(A)  \ 

0  ud^iX)  -V'^iX) 

\  pc^d^{X)  1 


ud'^{X)  ) 


is  a  generalized  symbol  of  the  discrete  operator  L**  (4.3).  The  entries  of  L*’(A)  are  generalized  symbols  of 
discrete  scalar  operators  defined  as  responses  of  these  operators  on  the  exponent  function  A-^  .  For  exam¬ 
ple,  the  generalized  symbol  of  the  central  second-order  difference  approximation  to  the  first  derivative  is 
=  d%X)V,d%X)  =  ^  (A- i).  Explicitly, 


T^(/i 


_  9‘(e*")  / -pc^9‘(e*")/i-|-n9"(e*")/2 


(4.9) 


-  _  -pc^9‘(e*")/i-|-n9"(e*")/2 
P  - 


Recall,  that  'D^(A)  =  JF^(A)  —  (v?d^{X)d'^{X)  —  c^(3'^(A))^) ,  and  JF^(A)  is  a  generalized  symbol  of  the  desired 
discrete  full-potential  operator.  The  choice  of  the  exact  solution  (4.1)  as  a  Fourier  mode  guarantees  that 
is  a  second-order  accurate  approximation  to  Qexact-  The  solution  (4.6),  (4.9)  satisfies  the  discrete  equation 
(4.3)  but  not  the  discrete  boundary  conditions. 

The  general  solution  of  the  homogeneous  system  of  equations  (3.6)  is  a  linear  combination  of  character¬ 
istic  solutions  2.k  in  the  form  zi^(j)  =  v/.Xl 


■^hom 


(4.10) 


=  Ip]" 


horn 

PLm 


=  XI  =  X 


where  are  linearly  independent  solutions  of  =  0  corresponding  to  the  roots  Xt  of  the  characteristic 

equation 

(4.11)  det  L‘’(A)  =  ud^{X)J^'^{X)  =  0. 

The  general  solution  of  the  discrete  equations  (4.3)  is 


(4.12) 

=  qhom  +  qpar  =  ^Ck'ik  + 

(  U  P 

4) 

K  p'"  ) 

k 

y  p  / 

where  ct  are  chosen  to  satisfy  boundary  conditions.  For  second-order  accuracy,  qfjo^  must  be  0(/i^)-small. 
The  discretization  error  function  is  defined  as 

(4-13)  q'‘-qXcf 
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5.  Analysis  of  Boundary  Conditions.  Well-posedness  of  a  linear  discrete  problem  is  equivalent  to 
solvability  and  solution  uniqueness  of  the  linear  system  of  equations  obtained  by  substituting  the  general 
solution,  q^,  into  the  equations  corresponding  to  (physical  and  numerical-closure)  boundary  conditions. 
The  unknowns  in  the  linear  system  are  the  coefficients  Ck-  This  condition  immediately  implies  that  the 
total  number  of  boundary-condition  equations  should  be  the  same  as  the  number  of  linearly  independent 
characteristic  solutions,  zu-  Overspecified  boundary  conditions  seemingly  violate  this  condition  but  often 
provide  a  well-posed  discrete  problem.  However,  as  we  will  show  later  in  Section  6,  the  set  of  overspecified 
boundary  conditions  is  equivalent  to  another  set  that  exactly  satisfies  the  relation  between  the  number  of 
boundary  conditions  and  the  number  of  linearly  independent  characteristic  solutions. 

The  number  of  linearly  independent  characteristic  solutions  (and  therefore  the  number  of  free  param¬ 
eters  Ck)  is  defined  by  the  length  in  mesh  sizes  (number  of  nodes  minus  1)  of  the  determinant  operator 
stencil.  The  stencil  should  be  considered  as  a  “maximal-length  footprint”  before  any  possible  cancellations 
occur.  The  stencil  actually  provides  a  more  detailed  information:  the  number  of  nodes  left  (right)  of  the 
center  indicates  the  total  number  of  required  (physical  and  numerical-closure)  boundary  conditions  at  the 
left  (right)  boundary.  Further  specification,  including  a  separate  count  for  physical  and  numerical-closure 
equations,  can  be  obtained  from  analysis  of  the  set  of  linearly  independent  characteristic  solutions.  For 
I-stable  discretizations,  the  characteristic  solutions  corresponding  to  |Ai;|  =  1  relate  to  the  solutions  of  the 
differential  equations  and  are  associated  with  physical  boundary  conditions.  The  characteristic  solutions 
corresponding  to  |Ai;|  <  1  relate  to  the  discrete  solutions  decreasing  exponentially  fast  as  functions  of  the 
distance  in  mesh  sizes  from  the  left  boundary.  These  characteristic  solutions  are  associated  with  (define 
the  number  of)  numerical-closure  equations  at  the  left  (inflow)  boundary.  Analogously,  the  characteristic 
solutions  corresponding  to  |Ai;|  >  1  relate  to  the  discrete  solutions  fast  decreasing  as  functions  of  the  distance 
in  mesh  sizes  from  the  right  boundary  and  define  the  number  of  numerical-closure  equations  at  the  right 
(outflow)  boundary.  Furthermore,  analysis  of  the  boundary-condition-sensitivity  (BCS)  matrix  relating  the 
coefficients  Ck  with  the  boundary  conditions  gives  a  precise  indication  of  B-stability  and  accuracy  order  of 
the  corresponding  solution.  For  B-stability,  all  the  coefficients  of  BCSmatrix  must  be  bounded  as  h  tends 
to  zero.  The  accuracy  order  can  be  estimated  from  the  amplitude  of  q^o^- 

The  following  example  illustrates  the  BCSmatrix  analysis  in  application  to  a  particular  second-order 
accurate  discretization  of  a  system  of  two  simple  uncoupled  equations.  The  BCSmatrix  analysis  of  the  dis¬ 
crete  problem  reveals  that  for  achieving  second-order  accuracy,  some  of  the  numerical-closure  equations  must 
be  second-order  accurate.  This  fact  contradicts  a  common  belief  that  one  can  keep  the  target  accuracy  order 
by  approximating  numerical-closure  equations  with  one-order  lower  accuracy  than  the  interior  equations. 

The  target  differential  equations  are 


(5.1) 


dxv  = 
dxW  =  g^. 


Let  the  exact  solution  of  the  differential  problem  defined  on  the  interval  a;  €  [0, 1]  be 

(  v{x)  \  _  (  Cx\ 


(5.2) 
then 

(5.3) 


w(x) 


9^{x) 

9^{x) 


C„ 


C, 

c„. 


lae 


and  the  boundary  conditions  associated  with  (5.1)  are 


(5.4) 


^'(0) 

w(0)  =  Cu,. 


The  second-order  accurate  discrete  scheme  employed  for  (5.1)  on  a  uniform  grid  with  mesh  size  h  is  defined 
as 


(5.5) 


y  d’‘  J  y  Wj  J  y  §^{jh)  j 


where  and  Wj  {j  =  0,1, . . .  ,N,N  =  1/h)  are  discrete  functions  representing  the  continuous  solutions  v{x) 
and  w(x),  and 

(5.6)  =  i  , 

^  (^i+2  “  +  6vj  -  4uj*_i  +  Vj_2'j  ■ 

The  determinant  of  the  matrix  in  (5.5)  is  a  seven-point  discrete  operator  with  the  maximal-length  footprint 
stencil 

(5.7)  ^  [  1  -8  22  -24  9  0  0  ] 

centered  at  the  underlined  position.  Thus,  a  set  of  boundary  conditions  for  the  discrete  equations  (5.5)  must 
include  four  boundary  conditions  at  the  left  boundary  and  two  boundary  conditions  at  the  right  boundary. 
The  characteristic  equation  for  (5.5)  is 


(5.8) 


^  -  8^  -I-  22^  -  24i  -I-  9  -I-  OA  -I-  OA^  =  0. 


Two  zero  coefficients  before  the  highest  powers  of  A  imply  that  there  are  two  infinity  roots  of  the  characteristic 
polynomial.  Equation  (5.8)  has  total  six  roots  Ai,2  =  1,  A3, 4  =  and  As^e  =  00.  A  set  of  linearly  independent 
characteristic  solutions  zi^(j)  =  is  given  by 

(5.9)  vi  = 
corresponding  to  Ai,2  =  1, 

(5.10)  ~  ^  1  ^ 

corresponding  to  A3,4  =  |.  Characteristic  solutions  corresponding  to  A5  and  Xq  exhibit  nonzero  values  only 
near  the  right  boundary  and  are  zero  in  the  interior  and  at  the  left  boundary.  Two  characteristic  solutions 
corresponding  to  As^e  =  00  are 


(5.11) 


and  zg 


0 

SN,j 


^k,j 


0,  if  j, 
1,  if  k  =  j. 


The  characteristic  solutions,  Zk,  are  normalized  to  satisfy  max \zk{j)\  =  0(1).  Two  sets  of  boundary  condi- 
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tions  for  the  discrete  system  (5.5)  are  tested: 
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Set  (A): 


(5.12) 


(i) 

Vo 

=  C,, 

(ii) 

Wo 

—  J 

(iii)a 

=  9^{h), 

(iv) 

i  (w?  - 

= 

(v) 

i  (f  ^AT-l  “  2w(v_2 

+  iWjv-s) 

=  9H^-h), 

(vi) 

i  (lW%  -2w%_^  + 

\wn--^ 

=  9\i) 

Set  (B)  is  similar  to  Set  (A)  with  condition  (iii)a  replaced  with  a  central  second-order  accurate  approx¬ 
imation 


(5.13) 


(iii)b  =9^ih)- 


The  equations  (i)  and  (ii)  correspond  to  the  physical  boundary  conditions,  all  other  equations  are  numerical- 
closure  equations.  The  only  difference  between  sets  (A)  and  (B)  is  the  approximation  order  of  the  numerical- 
closure  equation  (hi)  at  j  =  1.  The  general  solution  can  be  represented  in  a  form  similar  to  (4.12) 


(5.14) 


where  lo  =  ah, 


w- 


—  ^  Ck^k  + 


(5.15) 


V  =  iaC„ 


k^l 


w  =  la  ^  J 


and  the  generalized  symbols  of  discrete  first  and  fourth  derivatives  are 
(5.16) 


anA)=  + 

^^"HA)=  ^(^-4i+6-4A  +  A2). 


The  exact  solution  (5.2)  has  been  chosen  so  that  the  particular  solution,  (h,  w))^e*“'^ ,  of  the  problem  (5.5)  is 

a  second-order  accurate  approximation  to  {v,w)'^ .  Thus,  to  keep  second-order  accuracy,  the  part  of  (5.14) 
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related  to  the  homogeneous  solution  (  ^  cu^k)  must  be  0{h^).  Since  max \'Lk{j)\  =  0{l),k  =  1,  •  •  • ,  6,  the 

k^i  i 

amplitude  of  the  homogeneous  solution  is  determined  by  the  coefficients  c  =  (ci,  C2, . . . ,  ce)^. 

The  coefficients  c  relate  to  the  boundary  conditions  through  matrix  B. 


(5.17) 


B  c  -I-  d  =  0, 


where  vector  d  =  (di ,  ^2,  •  •  • ,  de)^  is  the  vector  of  items  independent  of  Ck  obtained  from  substitution  of  the 
general  solution  (5.14)  to  a  particular  set  of  boundary  conditions.  For  set  (A), 


(  1 
0 


(5.18) 


B  = 


0 
1 

0  0 
0  0 
0  0 
0  0 


■3/1 


16 


3h 

1 

1 

'3^-^h 


0 

0 

0 

0 

J_ 

2h 

_2 

h 


0  \ 
0 
0 
0 
0 


2h 


J 
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and 


(5.19) 


/  C,-v  \ 

Cw  -w 

g'^ih)  - 

^2(1  _  /j)  _  - Jllf - g*(Ar-l)cc) 


The  main  indication  that  the  chosen  set  of  boundary  conditions  is  legitimate  is  the  fact  that  matrix  B  is 
invertible,  i.e.,  detB  ^  0,  and  the  discrete  problem  (5.5),  (5.12)  is  well  posed.  Thus,  the  well-posedness 
requirement  (2)  to  practical  sets  of  boundary  conditions  (Section  1)  is  satisfied.  From  (5.17), 


(5.20) 


c  =  -B-^d. 


Normalization  of  d  entries  facilitating  the  analysis  of  B-stability  leads  to 

(5.21)  d  =  D  d. 


where  entries  of  the  diagonal  matrix  D  are  chosen  so  that  the  boundary  data  appear  in  d  with  0(1) 
coefficients.  For  many  practical  boundary-condition  sets,  the  coefficients  of  the  quantities  of  interest  are 
already  0(1)  in  d,  and  therefore  the  vectors  d  and  d  are  identical  (D  is  the  unity  matrix),  however,  the 
overspecified  boundary  conditions  result  in  distinguished  d  and  d.  The  BCSmatrix  is  defined  as 

(5.22)  BCS  =  B-iD-b 


In  general,  a  discrete  problem  is  B-stable  in  any  norm  (according  to  Definition  of  B-stability  in  Sec¬ 
tion  1),  if  all  the  entries  of  the  BCSmatrix  remain  bounded  as  h  tends  to  zero.  Even  with  some  unbounded 
BCSmatrix  entries  (tending  to  infinity  as  h  tends  to  zero),  the  problem  can  still  be  B-stable  in  some  com¬ 
mon  integral  norms,  such  as  or  L\.  This  is  the  case  when  the  unbounded  entries  are  located  in  the 
rows  of  the  BCSmatrix  corresponding  to  non-physical  characteristic  solutions  with  separated  from  1. 
Furthermore,  given  the  accuracy  of  boundary-condition  approximations,  the  BCSmatrix  predicts  the  size  of 
the  coefficients  Ck  and,  hence,  the  approximation  accuracy  of  the  discrete  solution.  Actually,  for  B-stability 
and  accuracy  predictions,  one  does  not  need  to  compute  the  exact  inverse  B“^;  it  is  enough  to  estimate  the 
asymptotic  behavior  of  the  entries  of  the  BCSmatrix  as  functions  of  h. 

In  our  example  (set  (A)), 


(5.23) 


BCS 


/  0(1)  0  0(/i)  0  0  0  \ 

0  0(1)  0(1)  0(/i)  0  0 

0  0  0(1)  0{h)  0  0 

0  0  0(1)  0  0  0 

0  0  o{h)  0  0{h)  0 

y  0  0  o{h)  0  0{h)  0{h)  ) 


where  o(K)  is  an  exponentially  small  function  of  h  (e.g.,  0^3“'^j).  The  discrete  problem  (5.5),  (5.12)  is 
B-stable  because  the  BCSmatrix  (5.23)  does  not  contain  any  entry  growing  as  h  tends  to  zero.  The  global 
accuracy  mainly  depends  on  ci  and  ci  because  Ai,2  =  1,  and  all  other  \]-  are  well  separated  from  zero.  From 
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Fig.  5.1.  Discretization  errors  in  v  and  w. 


the  first  row  of  (5.23),  the  accuracy  of  computing  the  coefficient  ci  depends  on  the  approximation  of  the 
boundary  conditions  (i)  and  (hi) .  If  boundary  condition  (i)  is  at  least  second-order  accurate  and  boundary 
condition  (iii)  is  at  least  first-order  accurate  (that  is  the  case  in  both  sets  (A)  and  (B))  then  the  coefficient 
Cl  is  computed  with  at  least  second-order  accuracy.  Similarly,  the  accuracy  of  computing  the  coefficient  C2 
depends  on  the  approximation  of  the  boundary  conditions  (ii),  (iii),  and  (iv).  From  the  second  row  of  (5.23), 
one  can  conclude  that  the  accuracy  order  in  computing  C2  usually  cannot  be  better  than  the  approximation 
order  for  the  boundary  condition  (iii).  In  set  (A),  the  boundary  condition  (iii)a  is  a  first-order  accurate 
approximation  to  the  first  equation  of  (5.5).  Therefore,  the  third  element  in  d  is  0(/i)-small.  Generally, 
when  a  set  of  boundary  condition  is  changed,  the  BCSmatrix  should  be  recomputed.  However,  if  two 
sets  differ  only  in  approximations  to  the  same  equations,  the  BCSmatrix  remains  unchanged.  In  set  (B), 
the  boundary  condition  (iii)b  is  a  second-order  accurate  approximation  to  the  first  equation  of  (5.5),  and 
dz  =  0{h^).  Thus,  one  can  expect  that,  with  set  (A)  of  boundary  conditions,  function  w  is  approximated 
with  first-order  accuracy,  and,  with  set  (B),  w  is  approximated  with  second-order  accuracy.  Function  v  is 
approximated  with  second-order  accuracy  in  any  case.  Figure  5.1  confirms  this  prediction.  On  the  figure, 
the  Loo-norm  of  discretization  errors  computed  for  =  2,Q!  =  Ttt  on  a  sequence  of  grids  with 

h  =  2~‘^,  2“®,  •  •  • ,  2“^®  is  shown. 

6.  Boundary  Conditions  for  Factorizable  Discretizations.  In  this  section,  the  boundary  condi¬ 
tions  for  two  versions  of  the  discrete  system  (4.3)  are  considered.  The  versions  differ  in  the  desired  form 
for  the  discrete  full-potential  operator  and,  therefore,  in  the  discrete  representation  for  the  operator  T>^. 
In  the  first  version,  the  discrete  full-potential  operator  employs  a  five-point  discrete  operator 
that  is  a  composition  of  second-order  accurate  upwind  and  downwind  difference  approximations  for  the  first 
derivative. 

h  +  4p^i  -  +  4p^i  -  \pU)  > 

#  {p'j+2  -  ^Pj+1  +  Qp'j  -  ^Pj-1  +  Pj-2)  ■ 

The  discrete  operator  used  in  the  second  version  is  a  three-point  central  approximation  as 

W  (Pj+i  - ‘^Pj  +  Pj-i)  ^ 

(p^+2  -  M+i  +  QPj  -  M-i  +  Pj-2)  ■ 


(6.2) 


gh  h  ^ 


(6.1) 


dLp’}  = 


XX^J 
rJl 
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All  other  operators  used  in  the  versions  of  (4.3)  are  the  same: 


dUp^h  ^ 

1_ 

h 

(l«i  - 

2<1 

+  \u>p_ 

2)  1 

d^p'p  = 

h 

-1)  ’ 

du^h  ^ 

h 

dV'i'  - 

- 

-2)  > 

d^u'p  = 

h 

-  in'* 
2 

-1)  ’ 

ddph  _ 

h 

{-M 

+  -  ip^+2)  , 

For  both  the  versions,  the  determinant  is  a  seven-point  upwind-biased  operator  implying  that  boundary 
conditions  must  contain  four  boundary  conditions  specified  at  the  inflow  boundary  and  two  boundary  con¬ 
ditions  specified  at  the  outflow  boundary.  Two  physical  boundary  conditions  (4.5)  are  defined  at  the  inflow 
boundary,  and  one  physical  boundary  condition  is  defined  at  the  outflow  boundary.  Thus,  two  inflow  and 
one  outflow  numerical-closure  equations  are  required.  In  examples  below,  the  following  sets  of  boundary 
conditions  are  analyzed: 

(A)  Overspecified  boundary  conditions,  where  values  oiu^i,UQ,u^ are  specified 
from  the  exact  continuous  solution  (4.1),  and  =  0.  The  total  number  of 

specified  values  is  twelve,  but  it  can  be  reduced  to  the  following  six  boundary  conditions: 


(i)  Mq  =Cu, 

(ii)  ud^u1  +  \d^p1  =fi(h), 

(iii)  ud^'tpi  —  =  0, 

(iv)  ud'^ip2  —  D'^p^  =  0, 

(v)  pc^d^u%_-^  +  +  ud‘^p%_i  =  /2(1  -  /i), 

(vi)  p%  =  Cpe*°^. 


The  equations  (ii) ,  (iii) ,  (iv) ,  and  (v) ,  used  in  the  boundary-condition  formulation  are  modified  from 
the  interior  equations  because  all  the  values  with  indexes  —1,  0,  N,  N  +  1  required  for  computation 
are  taken  from  the  overspecified  boundary  conditions  and  cannot  be  computed  (except  Uq  and  ppf) 
from  the  general  solution  representation  (4.12).  Note  also,  that  for  computation  of  (v)  one  should 
use  the  value  of  that,  in  turn,  computed  from  a  modified  equation  —  T>^p%_i  =  0. 

(B)  Practical  boundary  conditions.  Sets  of  practical  boundary  conditions  that  contain  only  physically 
available  boundary  data  are  described.  The  sets  differ  in  the  way  they  reconstruct  value  u'Pi. 


(i)  Uq  =Cu, 

(ii)  V'o  = 

(iii)  ul  (V:f  -  t/'o)  = 

(iv) b  7^?  (^2  -  3^5*  +  -  u(!.i)  =0, 

(v)  77?  {p%+i  -  2,p%  +  3p%_^  -  p%_2)  =  0, 

(vi)  p%  = 


The  conditions  (ii)  and  (iii)  imply  tplp  =  0.  One  can  interpret  the  boundary  conditions  (iv)b  and  (v) 
as  a  second-order  extrapolation  from  the  interior  for  the  value  of  u'Pi  and  p(v_|_i.  The  conditions  (iv)b 
and  (v)  are  numerical-closure  boundary  conditions  because  they  relate  outside  values  of  functions 
Uj  and  pj  with  inner  values.  These  two  conditions  are  equivalent  to  simultaneous  modification  of 
equations  ri(q*’)i  =  fi{h),r3{qfi)o  =  f2{0),ri{qfi)N  =  /i(l)  and  r3(q*’)Ar_i  =  /2(1  -  h).  (See  (3.6) 
for  definitions  of  ri;(q*’)j.) 
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(C)  Practical  boundary  conditions.  This  and  two  other  tested  sets  of  boundary  conditions  are  similar  to 
set  (B).  The  only  difference  is  the  condition  (iv)b  is  replaced  with 

(6.6)  (iv)c  {p^  -  +  3p'l;  -  p^)  =  0, 

The  condition  (iv)c  may  be  regarded  as  a  second-order  extrapolation  of  the  value  Pq  from  the  interior. 
The  interior  equation  r3(q*’)o  =  /2(0)  is  used  to  compute  The  condition  (iv)c  is  equivalent  to 
simultaneous  modification  of  equations  ri(q*’)i  =  /i(/i)  and  r2(q*')2  =  0. 

(D)  Practical  boundary  conditions. 

(6.7)  (iv)d  =«/i(0)-)/2(0), 

The  condition  (iv)d  reconstructs  from  a  second-order  accurate  central  discrete  approximation 
to  the  equation  {v?  —  c^)dxU  =  ufi  —  i/2  defined  at  the  inflow  boundary.  The  equation  is  derived 
by  manipulating  the  differential  equations  from  (2.5).  This  numerical-closure  boundary  condition  is 
equivalent  to  simultaneous  modification  of  equations  ri(q*’)i  =  /i(/i)  and  r3(q^)o  =  /2(0). 

(E)  Practical  boundary  conditions. 


(6.8) 


(iv)e  I  (mi  -  Mo)  +  ^  (iP2  -  iPo)  =hih), 


The  condition  (iv)e  is  shown  to  illustrate  the  importance  of  second-order  accuracy  in  approximating 
the  first  and  the  third  equations  in  (4.3)  at  the  inflow  boundary.  Simultaneously  with  (iv)e,  the 
interior  second-order  accurate  equation  ri(q*’)i  =  /i(/i)  is  used  to  recover 

6.1.  Factorizable  Scheme  #1:  5-point  Full-Potential  Operator.  Let  and  be  second-order 
upwind  and  downwind  discrete  operators,  respectively,  is  the  second-order  central  discretization,  and 
=  {u^  —  (?)d^d^.  The  generalized  symbols  are 


(6.9) 


H-f +2A-iA2), 

^'(A)=  i(-ix  +  iA), 

T'^{\)=  -c^)d^{\)d^{\), 

V'^(X)  =  c2  ((^"(A))^  -  a“(A)a'^(A))  =  -  4i  6  -  4A  -b  A^) . 


A  particular  solution  of  the  nonhomogeneous  problem  is  found  by  substitution  of  the  generalized  symbols  to 
(4.9).  The  characteristic  equation. 


(6.10)  (a“(A))^a'^(A)  =  0, 

has  six  roots  Ai,2,3  =  1,  A4,5  =  and  Xq  =  3. 

Before  presenting  the  general  solution,  let  us  introduce  some  auxiliary  symbol-like  functions  that  are 
required  for  defining  linearly  independent  characteristic  solutions. 


(6.11) 


5nA)= 

HH  +  ^A), 


The  six  linearly  independent  characteristic  solutions,  zi^(j)  =  v/^(j)Xl,  are  given  by: 
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( “  ^ 

(  jh 

\ 

(6.12) 

Vl  = 

0 

i 

,V2  =  0 

,  and  V3  = 

or}  (|J  -  1 

) 

1  0  J 

VW 

V  -jhpu 

/ 

corresponding  to  Ai,2,3  =  1, 


/ 

h  \  1 

(  h 

\ 

V4  = 

-pc^hd^{\i)  = 

|pc2 

\ 

0  )  ' 

V  0 

/ 

(6.13)  and 


(  jh 

] 

( 

jh 

V5  = 

P^^h  (fi  _  dc{\,)  -  jd^{\,)) 

= 

pc^  ( 

[-^^-l+jl) 

K  f^P^dHxl) 

) 

\ 

-jhpu 

/ 

corresponding  to  A4,5  =  and 


(6.14) 


V6  —  (Ae)  ^ 


h 


—hpu 


9"(A6) 

9"(A6) 


\ 

/  h  \ 

=  (As)-^ 

/ 

4  2  1 

-3PC 

V  -hlpu  J 

corresponding  to  Xq  =  3.  The  characteristic  solutions  are  normalized  to  satisfy  max|zi;(j)|  =  0(1)  for  h 

3 

tending  to  zero. 

For  overspecified  boundary  conditions  (set  (A)), 


(6.15)  d 


(  \ 

+  5pti) 

—  j.  3u^-4c^  x„h 

2h  "“at  ^  euh  °Pn+i 

V  ^p%  / 


/  1  0  0  0  0  0  \ 

0  /i  0  0  0  0 

0  0  /i^  0  0  0 

0  0  0  /i2  0  0 

0  0  0  0  /i  0 

^  0  0  0  0  0  1  y 


(6.16) 


and 


(6.17) 


/  Su^  \ 

—  ^(—4(5^0  +  <^V'-i)  +  c^(— 4(5po  +  ^P-i) 

—  ^SipQ  +  c^Spg 


^-Su% 

1  3u^-4c^ 
6u 

^P%+1 

\ 

) 

(  0(1) 

0(1) 

0(1) 

0(1) 

0(1) 

0{h) 

0(1) 

0(1) 

0(1) 

0(1) 

0(1) 

0{h) 

0(1) 

0(1) 

0(1) 

0(1) 

0(1) 

0(1) 

0(l/h) 

0(l/h) 

0(l/h) 

0(l/h) 

0(l/h) 

0(1) 

0(l/h) 

0(l/h) 

0(l/h) 

0(l/h) 

0(l/h) 

0(1) 

V  0(i/h) 

0(l/h) 

0(l/h) 

0(l/h) 

0(l/h) 

0(l/h) 
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where  Sqj  =  (q*'exact)i  “  (Q^par)ij  =  {Suj,S'tpj,Spj)'^.  The  solution  of  the  discrete  problem  (4.3)  with 
the  overspecified  boundary  conditions  is  obviously  B-unstable  according  to  Definition  of  B-stability  given  in 
the  introductory  Section  1.  Indeed,  upon  any  fixed  perturbation  of  the  overspecified  boundary  conditions,  the 
coefficients  C4,C5,  and  cq  grow  with  h  tending  to  zero.  However,  there  are  two  factors  that  allow  considering 
this  discrete  scheme  as  satisfactory: 

1)  The  unbounded  entries  of  the  BCSmatrix  affect  only  the  coefficients  of  the  “non-physical”  character¬ 
istic  solutions,  (i.e.,  the  characteristic  solutions  corresponding  to  Xt  with  amplitudes  separated  from 
unity)  and,  therefore,  decrease  exponentially  fast  away  from  boundaries  as  functions  of  the  mesh 
index.  So  the  discrete  solution  remains  bounded  in  any  common  integral  norm  (e.g.,  L2-norm). 

2)  Within  the  “non-physical”  characteristic  solutions,  only  the  amplitude  of  the  auxiliary  function 

is  unbounded  because  the  entries  of  V4,V5,  and  vg  related  to  physical  variables,  and  are 
0(/i)-small  in  comparison  with  the  entries  related  to 
Recall  that  the  values  of  are  0(/i^)-small.  Accordingly,  if  the  overspecified  boundary  conditions 
are  copied  precisely  from  the  exact  continuous  solution  (4.1),  the  discrete  solution  of  (4.3)  is  second-order 
accurate  for  the  physical  variables,  andp^,  in  any  norm;  the  Loo-norm  of  converges  to  zero  as  0{h),  and 
tp'^  is  0(/i^)-small  in  any  common  integral  norm.  We  found  experimentally,  that  with  exact  overspecification, 
only  Cq  behaves  as  0(h)-,  coefficients  C4  and  C5  are  still  0{hP).  This  better-than-expected  behavior  for  C4 
and  C5  is  explained  by  cancellations  in  computation  of  B“^d.  Upon  an  arbitrary  0{hP)  perturbation  of  the 
overspecified  data,  all  the  three  coefficients,  C4,C5,  and  cq,  become  0(h). 

Vectors  d  are  very  similar  for  the  sets  (B)  through  (E);  the  only  difference  is  the  fourth  element,  d^.  For 
set  (B), 

\ 


y 


(6.18) 


d  = 


/  d.  \ 
d2 
dz 

df^ 

V  ) 


Cu-'^ 

-4’ 


"-1 


-p- 


pi(JV  +  l)a,_3giJV 


utp- 

P 


-\-36 


i(N-l)u;  _  i(N-2)<. 


~h^ 

(f7p-p)e*^“ 


where  u,tp,  and  p  are  given  in  (4.9)  with  generalized  symbols  defined  in  (6.9). 


(6.19) 


XC)  _ 

—  P  ■■ 

df)  = 

'(E)  _ 


d 


(m/i(0)  -  i/2 

/i(h)-(nh^+pi^ 


2h  ’ 


The  BCSmatrices  for  tested  practical  boundary  conditions  fell  in  two  categories:  for  sets  (B)  and  (C), 


/  0(1)  0 

0  0(1) 

0  0(1) 

0  0 

0  0 

y  0  0 


O(h^) 

O(h^) 

0(h) 

O(h^) 

0(h) 

0(h?) 

0(h) 

0(h?) 

0(h) 

0(h^) 

o(h) 

o(h) 

o(h)  0  \ 

O(h^)  0(1) 

o(h)  0 

o(h)  0 

o(h)  0 

0(/i2)  0  j 
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(6.20) 


BCS  = 


Fig.  6.1.  Scheme  #  1:  Loo  norms  of  discretization  errors  in  u,  if,  and  p. 


and  for  sets  (D)  and  (E), 


/ 

0(1) 

0{h) 

om 

0{h) 

o{h) 

0 

\ 

0 

0(1) 

0{h) 

0(1) 

0{h^) 

0(1) 

0 

0(1) 

0{h) 

0(1) 

o{h) 

0 

0 

0(1) 

0{h) 

0(1) 

o{h) 

0 

0 

0(1) 

0{h) 

0(1) 

o(h) 

0 

V 

0 

o{h) 

o(h) 

o(h) 

0{h?) 

0 

J 

All  the  sets  of  practical  boundary  conditions  provide  solutions  to  the  discrete  problem  (4.3)  that  are 
B-stable.  Elements  di,d2,d^,  and  de  are  0(h?)\  d^  is  0(1);  ^4’^^  and  are  0(1);  ^4'°^  is  0(h?)\  and  ^4'^^ 
is  0{h).  Thus,  for  sets  (B)  through  (D),  the  computation  accuracy  of  coefficients  cu  is  at  least  second  order, 
i.e.,  the  Scheme  #  1  is  second-order  accurate  for  all  the  (physical  and  auxiliary)  variables  in  any  norm.  Eor 
set  (E),  first-order  accuracy  in  computing  C2,  C3,  C4,  and  C5  is  expected.  The  impact  of  characteristic  solutions 
Z2  and  Z3  is  global  for  all  variables  because  jh  =  0(1)  for  j  =  0{N)-,  therefore,  the  discretization  errors  in 
and  are  0{h)  in  any  norm. 

The  results  of  numerical  tests  are  shown  in  Eigure  6.1.  The  tests  have  been  performed  for  the  non- 
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dimensional  constant  coefficients  corresponding  to  a  Mach  number  M  =  0.5: 


u  =  cM, 


where  7  =  1.4.  The  parameters  of  the  exact  continuous  solution  defined  on  the  interval  x  €  [0, 1]  have  been 
chosen  as  Cu  =  l,Cp  =  2,  and  a  =  Itt.  The  Loo  norms  of  discretization  errors  have  been  measured  for  a 
sequence  of  uniform  grids  with  h  =  2~‘^,  2“®, . . . ,  2“^®.  The  results  confirm  exactly  the  predictions  made 
by  the  BCSanalysis.  An  important  observation  is  that  the  choice  of  boundary  conditions  strongly  affects 
the  absolute  value  of  the  discretization  errors.  The  gain  from  optimizing  the  set  of  boundary  conditions 
may  exceed  an  order  of  magnitude  in  solution  accuracy  on  a  given  grid.  Among  the  practical  boundary 
conditions,  set  (C)  seems  to  exhibit  the  minimum  discretization  errors  for  all  the  variables.  For  physical 
variables,  the  discretization  errors  with  overspecified  boundary  conditions  (set  (A))  is  minimal. 

6.2.  Factorizable  Scheme  3-point  Full-Potential  Operator.  Let  and  remain  the  same 
as  in  Factorizable  Scheme  #  1,  and  =  {u^  —  where  is  a  three-point  central  second-order 

accurate  approximation  to  the  second  derivative.  The  full-potential  factor  has  formally  a  three-point  stencil, 
however,  its  stencil  should  be  considered  as  five-point  long  with  zero  coefficients  (due  to  cancellations)  at 
the  outermost  nodes.  Therefore, 

(6.23)  =  1  _2  +  A  +  0A2j 

and 

(6.24)  V'^iX)  =  T\X)  -  [u‘^d%X)d\X)  -  c\d^{X)f)  =  (^1  _  4^  +  6  -  4A  +  j  . 

The  characteristic  equation 

(6.25)  a“(A)L''‘(A)  =  0 

has  six  roots  Ai^^s  =  1,  A4  =  |,  A5  =  0,  and  Xq  =  00.  The  same  five  sets  of  boundary  conditions  (sets  (A) 
through  (E))  are  analyzed  for  Factorizable  Scheme  #  2.  The  scheme  is  characterized  by  presence  of  zero 
and  infinite  values  of  Xk  ■ 

The  solution  representation  as  a  linear  combination  of  the  functions  is  relevant  only  for 

finite,  nonzero  eigenvalues.  For  zero  eigenvalue,  the  corresponding  characteristic  solutions  are  localized  at 
the  inflow  boundary,  i.e.,  they  exhibit  nonzero  values  at  the  inflow  and  are  zero  in  the  interior  and  at  the 
outflow  boundary.  By  analogy,  characteristic  solutions  corresponding  to  infinite  eigenvalue  are  localized  at 
the  outflow  boundary,  i.e.,  they  are  nonzero  only  at  some  locations  in  the  vicinity  of  the  outflow.  Infinite 
values  of  Xk  have  already  appeared  in  the  example  in  Section  5.  However,  in  that  example,  the  characteristic 
solutions  related  to  infinite  Xk  affected  only  one  variable,  w^.  In  general,  each  characteristic  solution  may 
affect  all  the  solution  variables.  In  such  a  case,  the  nonzero  values  for  primitive  variables  corresponding  to 
a  localized  characteristic  solution  are  usually  featured  at  different  j-locations.  The  exact  placements  of  the 
nonzero  entries  are  dictated  by  the  three  requirements: 

(1)  The  corresponding  matrix  B  relating  the  coefficients  of  characteristic  solutions  to  the  boundary 
conditions  is  invertible. 
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(2)  The  solution  values  outside  of  the  limits  defined  by  the  placements  are  not  needed  for  computing 
residuals  of  the  interior  equations. 

(3)  Multiplication  of  the  localized  characteristic  solution  by  a  constant  does  not  change  residuals  in  the 
interior  equations. 

Four  linearly  independent  characteristic  solutions  corresponding  to  finite  (nonzero)  eigenvalues  can  be 
found  in  the  usual  form 


(6.26)  vi  = 

correspond  to  Ai,2,3  =  1,  and 

(6.27) 


(  \ 

,  and  V3  = 

p{u‘^-c^) 

V  -jhpu  J 

V4 


/  h 

—hpc^d*^{Xi) 
\  0 


/ 

V 


corresponds  to  A4  =  |. 

The  placement  of  nonzero  values  of  the  localized  characteristic  solutions  depends  on  the  choice  of  bound¬ 
ary  condition  set.  For  overspecified  boundary  conditions  (Set  (A)),  the  characteristic  solution,  Z5,  localized 
at  the  inflow  (i.e.,  corresponding  to  A5  =  0)  is 


(6.28) 


Z5(i) 


/  ^  \ 

\  puhS\j  J 


and  the  characteristic  solution,  zq,  localized  at  the  outflow  (i.e.,  corresponding  to  Xq  =  00)  is 


(6.29) 


Z6(i) 


(  hSN-i,j  ^ 

\  —3puhSN,j  J 


For  practical  boundary  conditions  (sets  (B)  through  (E)),  the  locations  of  nonzero  entries  are  moving  one 
mesh  size  outward,  i.e.. 


(6.30) 


and 


(6.31) 


Z5(i)  = 


^  2 


'^O.i 


\  puh  (5o,i  / 


Z6(i)  = 


h  5  Nj 

S^_ 
—3puhSN+i, 


i.i 


For  overspecified  boundary  conditions  (set  (A)), 


(6.32)  D  = 


f  1 
0 
0 
0 
0 
VO 


0 

h 

0 

0 

0 

0 


0 

0 

0 

0 

0 


0 

0 

0 

0 

0 


0 

0 

0 

0 

h 

0 


0  \ 
0 
0 
0 
0 
1/ 


/ 


d  =  - 


Su^ 


2ph^ 


+ 


~hr~ 


■{-ASpq  +  Sp'l-^) 

Spo 


2h  mh  °p 


Quh 
'N 


^iV+1 
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(6.33) 


and 


(6.34) 


f  \ 

-ISu'l,  + 

—  ^(—4(5^0  +  <^V'-i)  +  (3w^  +  c^)(— 4(5po  +  ^P-i) 

—  ^StpQ  +  (3u^  +  c^)Spq 

-P^Xo.h  _  9sf+4cf.  c  ft 
2  "“jv  eu  ^Pn+1 

V  ^Pn 


/  0(1) 

0(1) 

0(1) 

0(1) 

0(h) 

0(h)  \ 

0(1) 

0(1) 

0(1) 

0(1) 

0(h) 

0(h) 

0(1) 

0(1) 

0(1) 

0(1) 

0(1) 

0(1) 

0{l/h) 

0(l/h) 

0(l/h) 

0(l/h) 

0(1) 

0(1) 

0(l/h) 

0(l/h) 

0(l/h) 

0(l/h) 

0(1) 

0(1) 

V  0(i/h) 

0(l/h) 

0(l/h) 

0(l/h) 

0(l/h) 

o(i/h)  y 

The  B-stability  and  accuracy  estimations  derived  from  the  BCSanalysis  of  the  overspecified  boundary 
conditions  are  very  similar  to  those  obtained  for  Factorizable  Scheme  #  1: 

1)  The  scheme  is  B-unstable  in  a  strict  sense,  but  the  instability  affects  only  auxiliary  variable  and 
only  in  a  0(/i)-small  vicinity  of  the  boundary.  The  scheme  is  B-stable  in  common  integral  norms. 

2)  With  the  exactly  overspecified  boundary  conditions,  the  zone  of  0{h)  accuracy  for  is  two  points 
adjacent  to  the  outflow  boundary  {tp%_2  and  V'at-i)-  All  other  solution  values  converge  with  the 
second  order  as  h  tends  to  zero.  The  absence  of  a  first-order  convergence  zone  at  the  inflow  is  a 
result  of  cancellations  implying  0(/i^ )-small  coefficients  C4  and  C5. 

3)  Upon  an  arbitrary  0{hP)  perturbation  of  the  overspecified  data,  the  coefficients  C4  and  C5  become 
0{h). 

For  practical  boundary  conditions,  vectors  d  are  the  same  as  (6.18),  (6.19).  The  BCSmatrices, 


(6.35) 


BCS 


for  sets  (B)  and  (C),  and 


/  0(1)  0  0(/i2)  0{h^)  0  0  \ 

0  0(1)  0{h)  0(/i2)  0  0(1) 

0  0(1)  0{h)  0(/i2)  0  0 

0  0  0{h)  0(/i2)  0  0 

0  0  0{h)  0(/i2)  0  0 

(^  0  0  0  0  0(/i2)  0  y 


(6.36) 


/  0(1)  0  0(/i2)  0{h)  0  0  \ 

0  0(1)  0{h)  0(1)  0  0(1) 

0  0(1)  0{h)  0(1)  0  0 

0  0  0{h)  0(1)  0  0 

0  0  0{h)  0(1)  0  0 

(^  0  0  0  0  0(/i2)  0  y 


for  sets  (D)  and  (E),  guarantee  B-stability  for  discrete  problem  (4.3)  and  second-order  accuracy  for  sets  (B) 
through  (D).  With  set  (E),  the  solutions  converge  with  first-order  rate. 

Numerical  tests  have  been  performed  for  Scheme  #  2  with  the  same  parameters  (6.22)  and  the  same  exact 
solution  {Cu  =  l,Op  =  2,  and  a  =  Ttt)  as  for  Scheme  #  1  in  Section  6.1.  Convergence  of  the  discretization 
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m  =  a  h 


Fig.  6.2.  Scheme  #  2:  Leo  norms  of  discretization  errors  in  u,  if,  and  p. 


errors  in  the  Loo  norm  is  demonstrated  in  Figure  6.2.  The  qualitative  results  confirm  the  predictions  of  the 
BCSanalysis  and  are  basically  the  same  as  in  case  of  Scheme  #  1: 

1)  The  Loo  norm  of  discretization  errors  in  physical  variables,  and  is  0{h^)  for  sets  (A)  through 
(D)  of  boundary  conditions  and  is  0{h)  for  set  (E). 

2)  The  Loo  norm  of  discretization  errors  in  the  auxiliary  variable,  is  0{h^)  for  sets  (B),  (C),  and 
(D)  and  is  0{h)  for  sets  (A)  and  (E);  for  set  (A),  is  0(/i^)-small  in  any  common  integral  norm. 

3)  The  minimum  discretization  errors  for  physical  variables  are  achieved  with  the  overspecified  bound¬ 
ary  conditions. 

4)  Set  (C)  is  the  best  among  practical  boundary  conditions  and  very  competitive  with  the  overspecified 
boundary  conditions. 

Quantitative  comparisons  between  Schemes  #  1  and  #  2  reveal  that  the  discretization  errors  in  physical 
variables,  and  p^,  with  Scheme  #  2  are  about  two  times  smaller  on  the  same  grids. 

7.  Numerical  Tests.  This  section  reports  results  of  numerical  experiments  performed  for  quasi-one- 
dimensional  subsonic  flow  in  a  convergent-divergent  channel.  The  nonconservative  nonlinear  differential 
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equations  describing  the  problem  are  defined  as 


udxu  +  =  0, 

(7.1)  -fpdxU  +  udxP  =  --fpu^, 

(-/ -  l)edxU  +  udxe  = 


where  7  =  1.4  and  <7(2;) 


(7.2) 


1  —  0.8a;(l  —  x)  is  the  area  distribution  term.  The  physical  boundary  conditions 


u{0) 

P(l) 

e(0) 


M(0)  [- 

M- 

7  [i 


1+2 


+M2(0)- 


+M2(0)7|r 
1  +  2^ 


7(7- 


1 _ _ 

r-l)  [1 


+M2(0)- 


0.53452, 

1.13984, 

2.04082 


corresponding  to  a  fiow  with  constant  entropy  and  a  Mach  number  of  0.5  at  inflow  and  outflow,  i.e.,  M(0)  = 
M(l)  =  0.5;  the  equations  are  nondimensionalized  by  density  and  speed  of  sound  at  the  sonic  condition. 
The  Mach  number  distribution  in  the  exact  solution  is  shown  on  Figure  7.1. 


Exact  Solution 


Fig.  7.1.  The  Mach  number  distribution  in  the  exact  solution 


Discrete  schemes  approximating  (7.1)  in  the  interior  belong  to  the  family  of  factorizable  schemes  de¬ 
scribed  as 


(7.3) 


uhd^^ph  _  pyhph 
'YPjd'^Uj  +  4>j  +  Ujd^^Pj 

(7  —  l)ejd'^Uj  +  Ujd^ej 


0, 

0, 


h  h  ^ 

-^PjUj- 


-(7-  l)eX 


where  cj*  =  a{jh).  The  boundary  conditions  related  to  the  fourth  (energy)  equation  are  the  physical 
boundary  condition  Cq  =  e(0)  and  a  numerical-closure  equation  corresponding  to  the  first-order  upwind 
discretization  in  the  discrete  equation  defined  at  j  =  1.  Six  different  discrete  scheme  have  been  tested 
for  the  other  three  equations: 

1.  Factorizable  Scheme  #  2  with  a  three-point  discretization  for  the  full-potential  factor  and  overspec¬ 
ified  boundary  conditions. 
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2.  Factorizable  Scheme  #  1  with  a  five-point  discretization  for  the  full-potential  factor  and  overspecified 
boundary  conditions. 

3.  Factorizable  Scheme  #  2  with  a  three-point  discretization  for  the  full-potential  factor  and  set  (C) 
of  practical  boundary  conditions. 

4.  Factorizable  Scheme  #  1  with  a  five-point  discretization  for  the  full-potential  factor  and  set  (C)  of 
practical  boundary  conditions. 

5.  Factorizable  Scheme  #  2  with  a  three-point  discretization  for  the  full-potential  factor  and  set  (B) 
of  practical  boundary  conditions. 

6.  Factorizable  Scheme  #  1  with  a  five-point  discretization  for  the  full-potential  factor  and  set  (B)  of 
practical  boundary  conditions. 
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Fig.  7.2.  Nonlinear  problem:  Loo  norms  of  discretization  errors  in  u,  p,  and  e. 


The  Loo  norms  of  discretization  errors  are  shown  on  Figure  7.2.  In  these  tests  performed  for  the 
nonlinear  nonconservative  Euler  equations,  the  discretization-error  history  demonstrated  by  Scheme  #  2  (3- 
point  approximation  for  the  full-potential  operator)  with  Set  (C)  of  practical  boundary  conditions  is  superior 
comparing  to  all  other  cases  of  practical  boundary  conditions  and  very  close  to  the  overspecified  boundary 
conditions  (Set  (A)). 
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8.  Conclusions  and  Further  Research.  Formulation  of  boundary  conditions  is  a  very  important 
step  in  design  of  discretization  schemes.  The  BCSanalysis  proved  to  be  a  very  accurate  and  reliable  tool 
in  predicting  the  effect  discrete  boundary  conditions  have  on  B-stability  and  accuracy  of  discrete  solutions. 
A  general  methodology  for  applying  the  BCSanalysis  to  discretized  well-posed  constant-coefficient  problem 
with  a  given  set  of  boundary  conditions  includes  the  following  steps: 

1.  Assume  the  exact  solution  of  the  target  differential  problem  to  be  a  Fourier  component  and  compute 
right-hand  side  and  boundary  data. 

2.  Compute  the  discrete  determinant  operator.  The  stencil  of  the  determinant  operator  determines  the 
required  number  of  boundary  conditions  at  the  inflow  and  outflow  boundaries. 

3.  If  the  total  number  of  boundary  conditions  (physical  boundary  conditions  and  numerical-closure 
equations)  exceeds  the  required  number,  an  equivalent  set  of  boundary  conditions  should  be  formu¬ 
lated. 

4.  Compute  the  roots  of  the  characteristic  polynomial. 

5.  Find  the  normalized  characteristic  solutions. 

6.  Find  particular  and  general  solutions  of  the  discrete  problem. 

7.  Substitute  the  general  solution  into  the  tested  set  of  boundary  conditions  and  compute  matrix  B 
and  vector  d. 

8.  Compute  the  BCSmatrix  and  vector  d.  To  evaluate  the  asymptotic  behavior  of  the  BCSmatrix  as 
a  function  of  mesh  size  h,  one  can  invert  matrices  B  numerically  on  several  grids  with  different  h 
and,  then,  compute  entry-by-entry  ratios  of  the  numerical  inverses. 

Practical  boundary  conditions  considered  in  this  paper  for  second-order  accurate  factorizable  discretiza¬ 
tions  are  easy  to  implement  and  result  in  B-stable  accurate  discrete  solutions.  The  BCSanalysis  and  nu¬ 
merical  tests  have  shown  that  for  second-order  accurate  solutions,  some  of  the  numerical-closure  equations 
must  be  second  order.  Optimization  of  practical  boundary  conditions  may  result  in  significant  gains  (of  more 
than  an  order  of  magnitude)  in  the  accuracy  of  discrete  solutions  on  a  given  grid.  Furthermore,  practical 
boundary  conditions  may  also  have  a  positive  effect  on  efficiency  of  multigrid  solvers.  It  is  expected  that 
implementation  of  practical  boundary  conditions  together  with  appropriate  local  procedures  accompanying 
distributed  relaxation  should  result  in  Laplace-like  multigrid  convergence  for  the  Euler  system  of  equations: 

(1)  TME  is  expected  with  multigrid  solvers  employing  V-cycles  rather  than  modifications  of  more  ex¬ 
pensive  W-  or  E-cycles  as  in  [9,  12,  18,  19,  20]. 

(2)  Assuming  very  efficient  solutions  of  convection  operators  (e.g.,  by  marching),  the  overall  conver¬ 
gence  rates  per  an  Euler-system  distributed  relaxation  sweep  is  expected  to  be  the  same  as  the 
corresponding  per-relaxation  rates  for  the  Laplace  operator,  e.g.,  0.38  per  sweep  of  lexicographic 
Gauss-Seidel  relaxation  within  a  V(l,l)  cycle  in  a  one-dimensional  subsonic  Euler  problem  with  the 
second-order  accurate  three-point  discretization  for  the  full-potential  factor.  In  a  similar  setting 
with  overspecified  boundary  conditions,  TME  solvers  with  distributed  relaxation  and  W  or  E  cycles 
usually  demonstrate  convergence  rates  around  0.5  per  relaxation. 

Preliminary  tests  confirm  the  expected  efficiency  in  solution  of  constant-coefficient  model  problems  corre¬ 
sponding  to  the  Euler  system  of  equations. 

The  BCSanalysis  can  be  regarded  as  a  generalization  to  systems  of  equations  of  the  half-space  discrete 
analysis  [7,  10]  that  takes  both  the  (inflow  and  outflow)  boundaries  into  account.  This  connection  implies 
that  the  BCSanalysis  can  be  used  far  beyond  the  task  of  boundary  condition  evaluation,  e.g.,  for  analyzing 
iterative  solvers.  It  can  be  extended  to  multidimensional  problems  by  considering  the  target  discretization 
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on  a  layer  bounded  by  two  parallel  hyper-planes.  Solutions  in  the  hyper-planes  parallel  to  the  boundaries 
are  represented  by  one  Fourier  component  in  a  time.  In  this  way,  the  target  multidimensional  problem  is 
translated  into  a  one-dimensional  problem,  where  the  frequencies  of  the  Fourier  component  are  considered 
as  parameters. 

The  applications  of  the  BCSmatrix  analysis  are  extended  far  beyond  the  examples  reported  in  this  paper. 
It  is  a  very  instrumental  in  gaining  insights  about  global  effects  some  local  conditions  may  have  on  solutions. 
The  BCSmatrix  analysis  plays  a  central  role  in  adjusting  the  interior  boundary  conditions  for  the  local 
procedures  supplementing  distributed  relaxation.  The  main  requirements  for  these  boundary  conditions  is  to 
prevent  the  local  procedure  from  error  magnifications  because  of  the  erroneous  data  unavoidable  introduced 
at  the  interior  boundary.  The  relative  simplicity  and  effectiveness  of  the  BCSanalysis  promise  many  further 
important  applications. 
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